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ABSTRACT 

The hot Jupiter HD 209458b was observed during primary transit at 3.6, 4.5, 5.8 
and 8.0 /Ltm using the Infrared Array Camera (IRAC) on the Spitzer Space Telescope. 
We describe the procedures we adopted to correct for the systematic effects present in 
the IRAC data and the subsequent analysis. The lightcurves were fitted including limb 
darkening effects and fitted using Markov Chain Monte Carlo and prayer-bead Monte 
Carlo techniques, obtaining almost identical results. The final depth measurements 
obtained by a combined Markov Chain Monte Carlo fit are at 3.6 /im, 1.469 ±0.013 % 
and 1.448±0.013%; at 4.5 /xm, 1.478 ±0.017% ; at 5.8 fxm, 1.549±0.015% and at 8.0 
/zm 1.535±0.011 %. Our results clearly indicate the presence of water in the planetary 
atmosphere. Our broad band photometric measurements with IRAC prevent us from 
determining the additional presence of other other molecules such as CO, CO2 and 
methane for which spectroscopy is needed. While water vapour with a mixing ratio of 
10 -4 — 10~ 3 combined with thermal profiles retrieved from the day-side may provide 
a very good fit to our observations, this data set alone is unable to resolve completely 
the degeneracy between water abundance and atmospheric thermal profile. 

Key words: techniques: photometric — planets and satellites: general — planetary 
systems — occultations 



1 INTRODUCTION 

More than 420 exoplanets, i.e. planets orbiting a star other 
than our Sun, are now known thanks to indirect detection 
techniques (Schneider, 2009). In the first decade after the 
initial discovery of a hot Jupiter orbiting a solar like star in 
1995 (Mayor and Queloz, 1995), the task was to find more 
and more of these astronomical bodies. In recent years, at- 
tention has switched from finding planets to characterising 



them. Among the variety of exoplanets discovered, particu- 
lar attention is being devoted to those planets that transit 
their parent star, and whose presence can therefore be de- 
tected by a reduction in the brightness of the central star 
as the planet passes in front of it. Sixty- nine of the 420+ 
currently identified exoplanets are transiting planets, and 
for these objects planetary and orbital parameters such as 
radius, eccentricity, inclination, mass (given by radial veloc- 



2 Beaulieu, Kipping, Batista, Tinetti, Ribas et al, 



ity combined measurements) are known, allowing first order 
characterisation on the bulk composition and temperature. 
In particular, it is possible to exploit the wavelength de- 
pendence of this extinction to identify key chemical com- 
ponents in the planetary atmosphere (Seager and Sasselov, 
2000; Brown, 2001), which permits enormous possibilities 
for exoplanet characterisation. 

The extrasolar planet HD 209458b orbits a main se- 
quence G star at 0.046 AU (period 3.52 days). It was the 
first exoplanet for which repeated transits across the stel- 
lar disk were observed (~ 1.5% absorption; Charbonneau 
et al., 2000). Using radial velocity measurements (Mazeh et 
al., 2000), the planet's mass and radius were able to be de- 
termined (M p ~ 0.69 Mj up , R p ~ 1.4 Rjup), confirming the 
planet is a gas giant with one of the lowest densities so far 
discovered. Consequently it must possess a highly extended 
atmosphere making it one of the optimum candidates for 
observation using primary transit techniques, and it was in- 
deed the first exo-atmosphere probed successfully using this 
method in the visible (Charbonneau et al., 2002) and then 
in the infrared (Richardson et al., 2006). 

Following the work on HD 189733b, where the first de- 
tections of water vapour (Tinetti et al., 2007b; Beaulieu et 
al., 2008) and methane (Swain, Vasisht & Tinetti 2008) have 
been achieved, we were awarded 20 hours Director's Discre- 
tionary Time on Spitzer (PI Tinetti, WETWORLD, PID 
461) to probe the atmosphere of HD 209458b in primary 
transit in the four IRAC bands at 3.6, 4.5, 5.8 and 8 /im 
(channels 1 to 4 respectively). Water vapour was proposed 
to be present in the atmosphere of HD 209458b by Barman 
(2007), to fit the data recorded by Hubble-STIS in the vis- 
ible (Knutson et al., 2007). Also, water vapour combined 
with a thermal profile increasing with altitude was a reason- 
able explanation to fit the secondary transit photometric 
data observed in the mid-IR (Deming et al., 2005; Knut- 
son et al., 2007; Burrows et al, 2007). Our understanding of 
the thermal profile and composition has improved thanks to 
more recent secondary transit spectroscopic data in the near 
and mid-IR, indicative of the additional presence of methane 
and carbon dioxide in the atmosphere of HD 209458b (Swain 
et al., 2009b), confirmed by Madhusudhan N. & Seager S. 
(2010). 

Transmission and emission spectra probe different re- 
gions of a hot- Jupiter atmosphere, both longitudinally and 
vertically (Tinetti & Beaulieu, 2008). In particular, the mid- 
infrared primary transit observations described here allow us 
to probe the terminator region of HD 209458b between the 
bar and millibar level. 



2 OBSERVATIONS AND DATA REDUCTION 
2.1 Planning the observations 

Three HD 209458 primary transits were observed with the 
IRAC camera on board the Spitzer Space Telescope. Chan- 
nels 1 and 3 (3.6 and 5.8 fim) were observed at two epochs, 
on December 30, 2007 and July 18, 2008, and data were ob- 
tained using channels 2 and 4 (4.5 and 8 /im) on July 20, 
2008. Since HD 209458 is a GOV star with a 2MASS Ks 
magnitude of 6.3, the IRAC predicted fluxes are 878, 556, 
351 and 189 mjy in channels 1-4, respectively. For our ob- 
servations we required extremely high signal-to-noise as the 



modelled contribution to the absorption due to H2O was 
predicted to be a few times 10 -4 of the stellar flux. 

As with other Spitzer observations of transiting plan- 
ets, it is necessary to observe the target continuously with- 
out dithering, in order to be able to quantify optimally the 
systematic effects detailed below. 

- Flat-fielding errors are an important issue; observa- 
tions at different positions on the array effectuate systematic 
scatter in the photometric data that can potentially swamp 
the signal that we are looking for. 

- The amount of light detected in channels 1 and 2 
shows variability that depends on the relative position of 
the source with respect to the pixel centre (labelled the pixel 
phase effect). The time scale of this variation is of the or- 
der of 50 minutes. These effects are well known and docu- 
mented in the IRAC Data Handbook and also discussed by 
Morales-Calderon (2006), Beaulieu et al., (2008), Knutson 
et al., (2008), Agol et al., (2008). To first order, these are 
able to be corrected for using the prescription of Morales- 
Calderon (2006). However, as the effects are variable across 
the array, ultimately they have to be estimated from the 
data themselves. 

- In channels 3 and 4 there are only minor pixel phase 
effects, but a variation of the response of the pixels to a long 
period of illumination and latent build-up effect impinge on 
the 5.8 and 8.0 /tm observations, respectively. 

- Wc obtained a slightly longer 'pre-transit' data set, 
in order to allow the satellite settle in a 'repeatable' jitter 
pattern and a shorter post-transit data set. The time scale 
of the pixel-phase effect being of the order of 50 minutes, 
we chose 120 min of pre- and 80 min of post-transit data 
baseline. 

It is important to note that the ~ 184 minute transit of 
HD 209458b means that our data contain three full cycles 
of the pixel phase variation in the transit itself, giving an 
excellent opportunity to have a full control on the behaviour 
of the systematic effects by evaluating them both in and 
outside the transit. 

Our observations employed the IRAC 0.4/2 second stel- 
lar photometry mode. Using the regular Astronomical Ob- 
servation Templates (AOTs), a total of two transits per field- 
of-view was required to achieve the desired sensitivity at 4.5 
and 5.8 /im (the arrays with the limiting sensitivity). Unfor- 
tunately, the AOTs as designed were not the most efficient 
way to perform this observation. Each stellar mode frame 
effectively incurs 8 seconds of overheads due to data trans- 
fer from the instrument to the spacecraft. As our observa- 
tions only required the data in the field-of-view with the 
star, it was possible to save both data volume by collecting 
data in only two channels and with a cadence of 4 seconds. 
Consequently, we designed a special engineering template 
(Instrument Engineering Request; IER) to optimise the ob- 
servations. IERs have been used successfully in other planet 
transit experiments (Charbonneau et al. 2005), and they can 
typically double the efficiency, and our IER enabled us to 
reduce the total required observing time for all four channels 
to only 13.4 hours. 

2.2 Data reduction and flux measurements 

We used the flat-fielded, cosmic-ray-corrected and flux cali- 
brated data files provided by the Spitzer pipeline. Each chan- 
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Figure 1. Raw photometric data for 3.6 fxm (epoch 1 and 2), 5.8 fim (epoch 1 and 2), 4-5 p.m and 8 fim obtained with IRAC. Each 
sub-panel has the same structure showing from top to bottom: the variation of the centroid position in X,in Y, and lastly the distance of 
the centroid from the lower left corner in the pixel (called the pixel phase, that can also be seen as the pointing error temporal amplitude). 
The lowest panel of each plot is the primary transit, and over-plotted the 50-point median-stack smoothing. They provide a synoptic 
view of the systematic trends present in IRAC primary transit data. 
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nel has been treated separately. We measured the flux of the 
target on each image using the version 2.5.0 of the SExtrac- 
tor package (Bertin & Arnouts, 1996), with a standard set of 
parameters for Spitzer (Infrared Array Camera Data Hand- 
book, 2006). The centroid determination was achieved with 
PSF fitting. We performed both aperture photometry, and 
PSF fitting photometry. In Fig.l, for all the six observed 
transits, we give the raw magnitude measurements (normal- 
ized using the post-transit observation), the variation of the 
centroid in X and Y axis, and the distance of the centroid 
from the lower left corner of the pixel (the pixel phase). A 
quick inspection shows that all observations contain corre- 
lated noise of different nature, as expected when using the 
IRAC camera. We discuss this phenomenon, and how we 
corrected for it, channel by channel, in the next section. 



3 ESTIMATION AND ATTENUATION OF 
CORRELATED NOISE 

3.1 Correcting the pixel phase effects 

It has been well-established that the IRAC channels exhibit 
pixel phase effects due to a combination of non-uniform re- 
sponse function within each pixel and very small pointing 
variations (Morales-Calderon et al., 2006, Beaulieu et al., 
2008). These effects are most prominent within the 3.6 /im 
and 4.5 /im photometry and to a lesser degree in the other 
two channels. We note that previous studies have not cor- 
rected for possible pixel phase effects at 5.8pm or 8pm, but 
in this work we evaluate the effectiveness of implementing 
it in all channels. 

Pixel-phase information is retrieved by using SExtrac- 
tor's PSF fitting to obtain estimates of the X and Y pixel- 
phase for each exposure. In contrast, the flux for each ex- 
posure is obtained through aperture photometry since this 
offers substantially larger signal-to-noise compared to the 
PSF flux estimates. 

A typical procedure is to directly correlate the X and 
Y phases to the out-of-transit fluxes to some kind of 4 or 5 
parameter fit (Morales-Calderon et al., 2006, Beaulieu et al., 
2008, Knutson et al. 2008) and in this work we will adopt 
a similar approach. We note that the PSF-fitted estimates 
of X and Y exhibit significant scatter at the same level as 
the amplitude of the periodic variations in each. This scatter 
is caused predominantly by photon-noise slightly distorting 
the PSF shape in a random manner and thus causing the 
fitting algorithm to deviate from the true value. The pixel 
phase effect is physically induced by the spacecraft motion 
and so we only wish to correlate to this property, as opposed 
to the random photon-noise induced scatter of X and Y. In 
order to do this, we fit a smooth function through the pixel- 
variations themselves before attempting to correlate to the 
out-of-transit flux. 

Our analysis of the X and Y phases reveals a domi- 
nant ~1 hour period sinusoidal-like variation in X and Y, 
characteristic of small elliptical motion in Spitzer's pointing, 
with a more complex time trend overlaid. For each channel, 
we apply a non-linear regression of a sinusoidal wave to the 
phases, in order to determine the best-fit dominant period, 
Pphase (typically close to one hour). We then calculate the 
median of the phases from the i th data point to the j th , 
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Figure 2. Zoom on the IRAC 5.8 fim observations to show the 
systematic trends that are present. First and second epoch, in the 
upper and lower panels respectively. We show the data and the 
50-point median-stack smoothing. For the second epoch notice 
the change of behaviour around 2800 seconds, indicated by the 
vertical line. Note that the behavior after 2800 sec is different 
between the two epochs. 



where tj — U + Pphase, (where tk is the time-stamp of the 
fc th exposure) and repeat from i = 1 up to the end of the 
data list. This moving- window- function essentially purges 
the dominant period from the phases and thus allows us 
to obtain a robust determination of the second-order phase 
variations, which may then be fitted for using a polynomial, 
of orders varying from 2 to 4 depending on the degree of 
curvature in the resultant phase trends. 

We have now calculated the function which describes 
the pixel phase variation of X and Y with respect to time, 
as induced by spacecraft motion. This function is then corre- 
lated to the actual out-of-transit photometry to find a fit to 
the function a + bX(t) + cY(t) + d[X{t)f + e[Y (t)] 2 . We find 
including an additional cross-term does not further improve 
the pixel-phase-effect attenuation. 

For 3.6pm (epochs 1 and 2) and 4.5pm, we removed 
pixel-phase effects of r.m.s. amplitude 0.49, 1.51 and 0.57 
mmag respectively, over the standard 8.4 second cadence. 
The second epoch at 3.6pm, is particularly polluted by pixel 
phase response, possibly due to a large inhomogeneity in re- 
sponse close to the PSF centroid position (pixel 131,128 of 
the detector). Repeating the process for the remaining chan- 
nels (after the other systematic effects were removed first, 
see next sections for details), we are able to remove 0.29 and 
0.24 mmag for 5.8pm (epoch 2) and 8pm respectively. Thus 
the pixel-phase induced variations are half of the minimum 
variations founds at 3.6pm and 4.5pm. 

3.2 Correcting systematic trends at 5.8 pm 

In the exoplanet community, at least two different meth- 
ods have been proposed to correct for the systematic ef- 
fects observed at 5.8pm, characterized by a large change 
in flux near the commencement of the observations. One 
frequently-adopted proceedure adopted is to discard the 
first ~30 min of observations (Knutson et al., 2008, Char- 
bonneau et al., 2008) and then de-trend the remaining 
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Table 1. Noise properties and effects of pixel- phase effect attenuation on each IRAC channel. 





3.6 ixm (epoch 1) 


3.6 fim (epoch 2) 


4.5 (im 


8.0 fim 


5.8 fim (epoch 2) 


5.8 fim (epoch 1) 


Before correction 


Baseline r .m.s. / Wag 


3.56005 


3.8848 


4.93071 


3.26671 


4.31389 


4.82004 


in % above photon noise 


75.5011 


92.3528 


119.734 


79.3238 


57.5393 


76.1831 


After correction 


Baseline r.m.s./mmag 


3.52621 


3.57796 


4.8974 


3.25821 


4.3042 


4.7886 


in % above photon noise 


73.8329 


77.1596 


118.249 


78.8573 


57.1856 


75.0341 


Noise removed/mmag 


0.489685 


1.51325 


0.572162 


0.235493 


0.288931 


0.549594 



data. For example, in the case of HD 189733b primary 
transit observations, Beaulieu et al. (2008) removed the 
first 20 min, and then applied a linear correction. Another 
method proposed by Desert et al. (2009) involves not ex- 
cluding these first 20 minutes but attempt to correct the 
data using a logarithmic parameterisation (see sec. 4.4): 
a + bt + clog(t — to) + d(log(t — to)) 2 . Employing different cor- 
rective procedures will undoubtedly yield significantly differ- 
ent transit parameters and so we must carefully consider the 
effect of each proposed correction. 

The most intuitive starting point is a visual inspection 
of the flux time series for our two measurements at 5.8(j,m. In 
figure[3] we exclude the transit event and show the behaviour 
of the out-of-transit flux (the baseline) with an overlaid 50- 
point median-smoothing as a visual guide. The first epoch 
exhibits a clear discontinuity between the photometry in the 
region ^ t < 5500 seconds and the subsequent data. The 
behaviour of this initial photometry does not match a 'lin- 
ear drift, a 'ramp' style-effect or any commonly employed 
analytic form. The origin of the observed behaviour is un- 
clear and is present in many different trial aperture sizes, 
between 2.5 to 20 pixels radius, suggesting an instrumental 
effect located either very close to the centroid position or 
globally across the detector array. 

Repeating the visual inspection for the second epoch, 
we observe a less pronounced version of this behaviour in 
the region ^ t < 2800 seconds (note that this behaviour is 
not seen in any other channels). However, the effect is osten- 
sibly sufficiently small that we cannot claim it is the same 
behaviour from a visual inspection of the time series alone. 
Therefore, we require a more in-depth analysis to provide a 
conclusion as to whether the systematic behaviours in epoch 
1 and epoch 2 are the same. In order to understand what 
kind of analysis this should be, we need to explicity qualify 
the question we are trying to answer. 

The difference between the truncation + linear trend 
versus the logarithmic correction can be summarized by one 
key point: the former proposes that the initial data is in- 
coherent with the latter data and cannot be characterized 
by a smooth analytic function. The latter works under the 
hypothesis that the entire time series is following one single 
smooth analytic description. We therefore wish to under- 
stand whether the properties of a smooth analytic function 
are consistent with the properties of the observed time se- 
ries. This is the critical question which we must answer. 

One key property of the smooth analytic, logarithmic 
function proposed by Desert et al. (2009), is that the dif- 




Time from first exposure/seconds 

Figure 3. Local gradient of each time stamp from the raw flux 
measurements obtained with IRAC at 5.8 fim for the two epochs 
of HD 209458b and HD 189733b. Note that the three exhibit 
similar behavior for the first 2000 sec. The first epoch for HD 
209458b has a larger amplitude of systematics, but the second 
epoch of HD 209458b and the observation of HD 189733b have 
remarkably similar behaviours. 



ferential of the function with respect time provides another 
smooth analytic function. In contrast, the truncation + lin- 
ear trend hypothesis postulates that since the initial data 
exhibits discontinuous behaviour, then the differential of this 
must also be discontinuous. So taking the differential of the 
time series will clearly resolve which hypothesis has the most 
supporting evidence. 

To achieve this goal, we first extract the uncorrected 
out-of-transit fluxes only and remove outliers for both epoch 
2 and epoch 1 using a median absolute deviation (MAD) 
analysis. We then create a moving 150-point window, in 
which we calculate the local gradient at each point. We 
do this by subtracting the median of the time stamps from 
each time stamp within a given window (to move the pivot 
along) and then performing a weighted linear regression. We 
define the weights as the square of the reciprocal of each 
flux measurement. In addition to the HD 209458b data pre- 
sented here, we perform the same process on the HD 189733b 
5.8/im data (used in Tinetti et al., 2007, Beaulieu et al., 
2008, Desert et al. 2009) for comparison giving us three data 
sets. Errors for each gradient stamp are computed using the 
weighted linear regression algorithm. 

In figure [3] we see all three local gradients plotted to- 
gether. Ostensibly, there seems to be strong correlations be- 
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Time from first exposure/seconds 

Figure 4. Local gradient of each time stamp from the raw flux 
measurements obtained with IRAC 5.8 fim (epoch 2) , computed 
using a linear regression of a moving 150-point window. Black 
indicates the observed local gradients, which differ greatly from 
those obtained using a logarithmic fit of the photometry (red). 
Notice the 1600 seconds peak/discontinuity. 



tween the three measurements, despite one of them being for 
a completely different star. In particular, there is a strong 
dip at around 2000 seconds after the first exposure in all 
three observations. In figure , we plot just the epoch 2 of 
HD 209458 and also overlay the local gradients obtained 
from a logarithmic fit of the baseline (equivalent to the first 
differential of this function with respect to time). It is clear 
that the logarithmic fit cannot explain the strong negative 
peak observed in the gradients data. Furthermore, the clear 
presence of discontinuous behaviour in the local gradients 
supports the hypothesis that no continuous analytic func- 
tion can correct this behaviour. 

Although the three measurements appear correlated, we 
may quantify these correlations. Comparing any two chan- 
nels, we define one as the reference data and one as the 
comparison data. We first ensure the minimum to maximum 
time stamps of both sets are the same by clipping the longer 
set appropriately. We then perform a linear interpolation of 
both the gradient measurements and the uncertaintities, for 
the comparison data. This allows us to accurately estimate 
the gradient values at like-for-like time stamps. Regenerat- 
ing the comparison gradients data using the interpolation 
function, we evaluate the correlation between the two using 
Pearson's correlation coefficient. We repeat the same process 
for randomly generated data with the same uncertainities 
and array length as the original. This is repeated 100,000 
times in order to estimate the expected correlations from 
random noise. 

Although they have been taken more than 7 months 
apart, we find epoch 1 and 2 for HD 209458 have corre- 
lation Corr(209458epoch2,209458epochl) = 0.678. The 10 5 
randomly generated noise values yield 0.001 ± 0.037. This 
makes the correlation significant at the 18.4-cr level. Re- 
peating the exercise for epoch 2 of HD 209458 and for the 
observations of HD 189733b (taken 3 years apart), we have 
Corr(209458epoch2, 189733) = 0.658. The randomly gener- 
ated noise yields 0.002 ± 0.043, making the observed cor- 
relation significant at the 15.2-ct level. In conclusion, the 
correlations in the local gradient plots are highly significant 
even for observations separated by years on different stars. 



We therefore conclude the observed behaviour must be in- 
strumental effects for 5.8/im detector array itself. 

The largest feature is that of the 'negative spike' at 
around 2000 seconds. After this, all 3 observations exhibit 
variations consistent with that of a singular constant value 
i.e. a linear fit. The reduced \ 2 of these three channels may 
be computed both for all data and for those data after the 
negative spike. We find the values always decrease by ex- 
cluding the negative spike; quantitatively we have respective 
changes of 1.18 -> 0.58 for HD 209458 epoch 2, 2.99 -> 1.24 
for HD 209458 epoch 1 and 2.66 -> 1.51 for HD 189733. 
Therefore we can see that the instrumental systematic ef- 
fects of 5.8/im can be split into two regimes, pre and post 
spike. The pre-spike data exhibits discontinuous behaviour 
compared to the latter data and cannot be characterized by 
a smooth continuous function. The post-data conforms to a 
linear fit. 

We therefore conclude that an analysis of the differen- 
tial of the time series supports the hypothesis that the 5.8/im 
correction should be to remove the discontinuous data before 
the gradient spike and then use a linear fit for the remainder. 
It would therefore seem that at 5.8/mi the detector requires 
a certain amount of time to settle into a stable regime, as in- 
dicated also in earlier studies (Beaulieu et al., 2008, Knutson 
et al., 2008, Charbonneau et al., 2008). 

Despite the evidence from this gradients analysis, we 
may conceive of several other possible tests to be certain 
that the logarithmic correction not favoured by the data. 
Using the lightcurve fitting code described in §4.3, we fitted 
two possible systematic correction lightcurve: 1) a trunca- 
tion of the first 2800 seconds, followed by a linear fit to the 
remaining baseline data (previous examples Knutson et al. 
2007; Harrington et al. 2007; Beaulieu et al. 2007) 2) a log- 
arithmic fit to all baseline data (previous example Desert et 
al. 2009). We select several properties to compare these two 
possible corrections: 

(i) Adopting a baseline between 2814 ^ t ^ 7543 seconds 
(i.e. after the discontinuous behaviour) and 19605 ^ t ^ 
23974 seconds, constituting 1082 data points, we compute 
the x 2 f° r both the linear and the logarithmic fit. Despite 
using two extra free parameters, the logarithmic produces a 
larger \ 2 = 1358.1 compared to a linear fit with \ 2 = 1338.7 
(flux uncertainties based on photon noise only). 

(ii) We may also compare the \ 2 °f t ne entire lightcurve 
fit (using the model described in §4.3). In this case, we must 
scale the x 2 f° r a fair comparison since the linear fit uses 
fewer points due to the truncation procedure. Comparing 
the reduced \ 2 between the two corrections we find lower 
values for the linear fit again- 1.045 vs 1.014 or 1.031 vs 
1.000, depending whether we additionally correct for pixel 
phascQ- 

(iii) We use the fitted transit duration, T, defined by 
Carter et al. (2009) which was shown to be highly robust and 
non-degenerate. T is expected to be independent of wave- 
length as the only possible parameter which could vary is R* 
which is not expected to exhibit significant changes between 
different wavelengths. We therefore refit the Brown et al. 
(2001) HST lightcurve of HD 209458b, taken in the visible, 

1 Reduced \ 2 values have been rescaled so that lowest value is 
equal to unity 
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with the same model used here. We find a transit duration of 
95251 i| seconds. In comparison, correcting the second epoch 
at 5.8 nm. with a linear fit yields Tu n = 9518^55 seconds and 
with a logarithmic fit Ti og — 9546lg5 seconds. 

(iv) In Fig 4., the local gradients, as taken in 150-point 
bins, is compared to that expected from the logarithmic fit 
of the data. There is a very strong discrepancy between the 
data and the model before 2000 seconds. 

Thus, we find employing a logarithmic fit, with two ad- 
ditional free parameters, cannot be shown to offer any kind 
of improvement over the linear fit. The gradient analysis pre- 
sented above shows that the logarithmic parametrisation is 
not adapted. Moreover, it is disfavored by A\ 2 = 20. Since 
every single test performed has supported the truncation 
and linear trend correction, this method will be adopted at 
the preferred corrective procedure in our subsequent analy- 
sis. 

3.3 Correcting the ramp at 8 /im 

The ramp effect at 8pm is well documented and so too is 
the methodology for correcting this phenomenon (Agol et 
al., 2008 and references therein). Unlike the 5.8//m data, 
there are no known discontinuities in the time series and 
thus the correction may be achieved using a smooth analytic 
function. We fit a time trend to the out-of-transit data of the 
form a + bt + clog(t — to) + d(log(t — to)) 2 where to is chosen 
to be 30 seconds before the observations begin to prevent 
the function exploding at t = 0. 



4 FITTING THE TRANSIT LIGHT CURVES 

Among the 6 transit light curves, we have four of high quality 
with well understood and corrected systematic effects : the 
first epoch at 3.6 /im, 4.5 /jm, the second epoch at 5.8 /im 
and 8 ^m. The second epoch at 3.6 fim and the first epoch 
at 5.8 /im will be treated separately. 

4.1 Limb darkening 

Accurate limb darkening coefficients were calculated for each 
of the four IRAC bands. We adopted the following stellar 
properties: T cff = 6100 K, log 5 = 4.38, and [Fe/H] = 0. We 
employed the Kurucz (2006) atmosphere model database 
providing intensities at 17 emergent angles, which we in- 
terpolated linearly at the adopted T e ff and log<? values. 
The passband-convolved intensities at each of the emergent 
angles were calculated following the procedure in Claret 
(2000). To compute the coefficients we considered the fol- 
lowing expression: 

7(1) = ^~l^Ck{l-V ), 

where I is the intensity, fj, is the cosine of the emergent angle, 
and Cfc are the coefficients. The final coefficients resulted 
from a least squares singular value decomposition fit to 11 
of the 17 available emergent angles. The reason to eliminate 
6 of the angles is avoiding excessive weight on the stellar limb 
by using a uniform sampling (10 /i values from 0.1 to 1, plus 
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Table 2. Limb darkening coefficients. 



channel 


cl 


c2 


c3 


c4 


(3.6 fim) 


0.2670569 


0.1396675 


-0.1900802 


0.064018 


(4.5 fim) 


0.3325055 


-0.1999922 


0.1858255 


-0.0703259 


(5.8 ftm) 


0.3269256 


-0.2715499 


0.2258883 


-0.0684003 


(8 nm) 


0.2800222 


-0.2278080 


0.1451840 


-0.0273881 



ii = 0.05), as suggested by Dfaz-Cordoves et al. (1995). The 
coefficients are given in Table [2] 



4.2 Markov Chain Monte Carlo fit to the data 

We adopt the physical model of a transit light curve through 
the expressions of Mandel & Agol (2002) and orbital eccen- 
tricity using the equations of Kipping (2008). We sampled 
the parameter space with Markov Chain Monte Carlo codes 
(Doran & Muller 2004) originally developed for microlensing 
(Dong et al., 2008; Batista et al., 2009) and adapted to fit 
transit data. We first made an independent fit for 3.6 fim 
(epoch 1 and 2), 4.5 /jm , 5.8 fim (epoch 2) and 8 [im. We 
adopted a fixed value of period to be P = 3.524749 days 
(Knutson et al. 2007). For each channel, 5 parameters are 
fitted, namely the out-of-transit baseline, the orbital incli- 
nation i, the ratio between the orbital semi- major axis and 
the stellar radius a/R*, the ratio of radii, k, and the mid- 
time transit t c . We also permit the orbital eccentricity e and 
the position of periastron cj to move in a restricted range, 
corresponding to the best-fit values derived by Winn et al. 
(2005) including their error-bars. The five other parameters 
are free. The error-bars of the data have been rescaled to 
make the \ 2 P er degree of freedom equal to unity. The re- 
sults are shown in Table [3] 

As some physical parameters should be the same for all 
bands, we made a simultaneous fit to the best observations, 
namely 3.6 /J.m (epoch 1), 4.5 /im , 5.8 fim (epoch 2) and 
8 /jm, in which four parameters are shared by all channels: 
P, e, i, ui and a/R*. Three other parameters, k, t c and the 
baseline, are fitted independently for each band and are al- 
lowed to move within the range obtained in the individual 
fits. We decided to fit separately 3.6 \xm (epoch 2), forcing 
the four shared parameters to be equal to the values derived 
from the best fit with the four other channels. The results 
are shown in Table 0] 



4.3 Prayer-bead Monte Carlo fit to the data 

We also fitted all the transit data with the code used by 
Fossey et al. (2009), incorporating the effects of non-linear 
limb darkening through the expressions of Mandel & Agol 
(2002) and orbital eccentricity using the equations of Kip- 
ping (2008). We fixed the orbital eccentricity, e, and position 
of periastron, w, to the best-fit values derived by Winn et al. 
(2005), and adjusted k, a/R*, b, and t c to find a minimum 
in x 2 - Although the parameters a/R* and 6 show a degree of 
covariance, Carter et al. (2008) have shown that the transit 
duration, T, ratio-of-radii, k, and mid-transit time, t c , are 
non-degenerate parameters; these parameters are also less 
affected by systematic errors in orbital eccentricity and thus 
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Figure 5. Final light curves, best fit model and residuals at 3.6 [im (epoch 1 and 2), 4-5 fim, 5.8 fim (epoch 2) and 8 fim. In the 
residuals subpannel we will overplot the 50-point median-stack smoothing of the residuals. 



Table 3. Markov Chain Monte Carlo fit to individual primary transits observed by IRAC. We list all the fitted parameters (see the text 
for the description), and in particular the ratio- of -radii, k = R p /R,, the orbital semi-major axis divided by the stellar radius, a/R*, 
the orbital inclination, i and the mid-transit time, t c . 





3.6 fim (epoch 1 ) 


3.6 fim (epoch 2) 


4.5 (im 


5.8 fim (epoch 2) 


8 fim 


i(deg) 


87.00 ±0.11 


86.67 ±0.15 


86.87 ±0.10 


86.84 ±0.11 


86.37 ±0.13 


A/R, 


8.89 ±0.06 


8.84 ±0.10 


8.91 ±0.05 


8.84 ±0.07 


8.49 ± 0.08 


k = R p /R, 


0.120835 ± 0.00054 


0.120387 ±0.00053 


0.1218 ±0.00072 


0.1244 ±0.00059 


0.1240 ±0.00046 


k 2 = (R p /R,) 2 


1.460 ± 0.013% 


1.449 ± 0.013% 


1.4835 ± 0.017% 


1.547 ±0.015% 


1.538 ± 0.011% 



can be taken to be more reliably constrained than a/R*, b, 
or the inclination, i. 

We use the genetic algorithm pikaia (see Metcalfe & 
Charbonneau 2003) to find an initial, approximate solution, 
which is used as the starting point for a x 2 - mm i m i sa tion 
using the downhill-simplex AMOEBA algorithm (Press et al. 
1992). The initial best-fit parameters from AMOEBA are ran- 



domly perturbed by up to 40% of their value and refitted in 
100 trials to check the robustness of the best-fit solution. 

To obtain the final parameter uncertainties, we em- 
ploy a 'prayer-bead' Monte Carlo simulation of the unbinned 
data, as used by Gillon et al. (2007). Here, the set of residu- 
als from the best-fit solution is shifted by one data point and 
added to the best-fit transit model to generate a new data 
set, with the residual at the end of the data series wrapping 
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Table 4. Markov Chain Monte Carlo fit to IRAC data. The first column shows a join fit to the best four band observations, namely 
3.6 fim (epoch 1) , 4-5 fim , 5.8 fim (epoch 2) and 8 fim. Then we impose the parameters i, A/R*, e, ui and fit the ratio of the radii k 
and the mid transit time t c for the second epoch at 3.6 fim and the first epoch at 5.8 fim . 





3.6 fim (epoch 1 ) + 4.5 fim 
+ 5.8 fim (epoch 2) + 8 /im 


(3.6 fim) epoch 2 


(5.8 fim) epoch 1 


i(deg) 
a/R t 


86.76 ±0.10 
8.77 ±0.07 






k = (R p /R,) 
(3.6 fim) 
(4.5 fim) 
(5.8 fim) 
(8 fim) 


0.121215 ± 0.00054 
0.121568 ± 0.00072 
0.1244 ±0.00059 
0.12390 ± 0.00046 


0.120343 ±0.00053 


0.1246 ± 0.00095 


k 2 = (R p /R,) 2 
(3.6 fim) 
(4.5 fim) 
(5.8 fim) 
(8.0 fim) 


1.469 ± 0.013% 
1.478 ± 0.017% 
1.549 ±0.015% 
1.535 ± 0.011% 


1.448 ± 0.013% 


1.552 ±0.032% 



1.01 



D 



S 0.99 



0.005 



-0.005 
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1 t I ' 1 
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-200 -100 100 

Time from mid transit (miri) 

Figure 6. The uncorrected unbinned and binned (30-points) data 
from first epoch at 5.8 fim and the underlined model computed 
for the second epoch (corrected for systematics) . The lower pannel 
shows the residuals of the binned data (unbinned data ommitted 
for clarity). The shaded area is marking the second half of the 
transit and the post transit observations used in the fit. 



around to the beginning. The new data set is refitted, and 
the process repeated until the set of residuals has been cy- 
cled through the entire data series. This procedure has the 
advantage of preserving the structure of any residual cor- 
related noise within the light curve in each simulation. For 
the unbinned data we then have typically 2500-3000 samples 
from which the parameter uncertainties may be estimated, 
which we take to be the values comprising 68.3% of the sam- 
ple about the median of each parameter distribution. The 
median and uncertainties are compared to the fitted value 
in each case, to check the robustness of the simulations and 
to assign upper and lower limits on the parameters. In all 
cases, we found the difference between the median and the 
best-fit parameter was insignificant. 

Table 5 lists the fitted depth, ratio of radii, k, transit 
duration, T, orbital inclination, i, and a/R* from this fitting 
procedure, for each of the transits. 



4.4 The case of the epoch 1 of 5.8 ftm 

The first epoch of 5.8fim requires special consideration due 
to the extremely pronounced nature of the systematic errors 
for this data set. In §3.2 we demonstrated that the corrective 
procedure most consistent with the observational evidence 
is to the truncate the initial data exhibitting discontinous 
behaviour and then perform a linear correction through the 
remaining baseline data. Since the systematic effect is so pro- 
nounced here, we took a conservative approach by assuming 
the systematic may persist up the moment of mid-transit. 
We therefore exclude all data before the mid-transit and 
adopt the physical parameters (i, a/R t ) derived from the 
global MCMC fit and reported in Table 4, and fitted the 
baseline, the ratio of radii k and the mid transit time t c . 

In Fig.[6]we compare epoch 1 data and the model fitted 
on the data from the mid-transit. Inspection of the residuals 
after the mid-transit indicates a good fit to the data. The 
data from the first half of the observations show the uncor- 
rected systematic trends at work. It is clear that they are of 
different nature from the one of epoch 2. We add the mea- 
sured ratio of radii and depth in the last column of Table 
4, and last row of table 5. The results we therefore obtain 
from second epoch is consistent with the first epoch. 

As a final check of the procedure, we decided to treat 
also the second epoch at 5.8 ftm the same way. We take the 
uncorrected data, exclude the first half of the data up to the 
mid-transit, and fit the light curve. We report the measured 
depth by this procedure to be 1.540 ±0.029%. It is perfectly 
compatible with our complete fits reported in Tables 3, 4 
and 5. 



4.5 Sanity check: Grid calculations 

As a check on our methodology we imposed the physical 
parameters derived from Knutson et al. (2008) in the ap- 
proximation of circular orbit, and fitted for the baselines, 
ratio of radii, k, and mid-time transit, t c , using a simple y 2 
minimisation. This gradient based method explores the local 
minimum around the physical solution found by Knutson et 
al., (2008). By comparing the results to the ones obtained 
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Table 5. Best-fit transit depths, ratio of radii k, duration T, orbital semi-major axis divided by the stellar radius a/R t , inclination i, 
mid transit time t c found using the prayer-bead Monte Carlo fit method described in §^._Z 



band//xm 



(Rv/R-f 



Rp/ R* 



T/seconds a/R t 



3.6 /im (epoch 1) 

4.5 (im 
8.0 fim 

3.6 /im(epoch 2) 
5.8 /im(epoch 2) 



1 462+ ' 011 

1 449+0.010 

,+0.0099 
-0.0096 



1.542" 1 



-.+0.00045 
^-0. 00048 
.+0.00056 
± -0. 00056 
12403+ 00043 

U.12U38_ 00043 
o + 0.()(l()4() 
-0.00039 



0.12089] 
0.12174^ 



0.12416 1 



9581 
9437 
9580 
9358 
9517 



+57 
47 
+60 
51 
+58 
49 
+57 
48 
50 
54 



8.88 ±0.02 
9.09 ±0.01 
8.50 ±0.02 
8.86 ±0.02 
9.13 ±0.01 



86.99 ±0.18 
87.15 ±0.10 
86.32 ±0.20 
86.70 ±0.13 
87.22 ±0.12 



with the Monte Carlo methods, we found that the transit 
depths agree well within the error bars. By contrast the other 
parameters (inclination, a/R*) are highly degenerate. 

4.6 Influence of spots 

An effect to consider when comparing transit depth at dif- 
ferent wavelengths is the influence of stellar surface inhomo- 
geneities, i.e, star spots. Depending on the spot distribution, 
the occulted stellar area during the transit can be brighter 
or dimmer than the average photosphere. In the case of HD 
189733, a moderately active star with visible photometric 
variations of ~3 % (peak to peak), the differential effect in 
the IRAC 3.6 and 5.8 fim bands was evaluated by Beaulieu 
et al. (2008) to be below 0.01%. HD 209458 is a chromo- 
spherically inactive star with an estimated age close to that 
of the Sun (e.g., Mazeh et al. 2000; Cody & Sasselov 2002; 
Torres et al. 2008). It is thus reasonable to assume a level 
of photometric variations similar to that of the Sun, i.e., 
0.2-0.3% peak to peak (Frohlich & Lean 2004). Scaling the 
calculations carried out for HD 189733, the expected differ- 
ential effect of star spots on the IRAC bands for HD 209458 
is likely to be 10-20 times smaller, and therefore well below 
0.001%. Our calculations show that star spots have negligi- 
ble influence when compared with our measurement uncer- 
tainties (~0.011-0.017%); see Tables I3l4l and l5l 

4.7 Comments about different epochs at 3.6 and 
5.8 /im 

We asked for two epochs for HD 209458b at 3.6 and 5.8 /im 
with the prime intention of demonstrating the possibility of 
co-adding multiple epoch observations, and/or to be able 
to check for the variability in the system. The two epochs 
are separated by 7 months, and the observing setups are 
identical. 

Firstly, at 3.6 /jm the data are affected by systematic 
trends of the same nature due to the pixel scale effect. We 
notice a factor 3 in the amplitude of the systematic trends 
between the two epochs. We measure the two transit depth 
to be 1.469 ± 0.013% and 1.448 ± 0.013% respectively. The 
results are compatible between the two channels. 

Secondly, at 5.8 /jm the situation is more complex. The 
second epoch showed the expected behaviour, and we have 
been able to correct for systematics, and to fit it. For the 
first epoch, we choosed to discard the first half of the data, 
and fit the uncorrected remaining data set. We measure the 
two transit depth to be 1.552 ± 0.032% and 1.549 ± 0.015% 
respectively. The results are compatible between the two 
channels. 



Even when centering on the same pixels of the detector, 
observing the same target several months apart, different 
systematics are at work. It is clear that the different data 
sets should be analysed for systematics and then corrected 
individually. Then, multiple epoch can be compared and/or 
added. 

4.8 Results 

We have chosen three approaches to fit the data, i.e., grid 
calculations with ephemeris from Knutson et al., 2008, 
Markov Chain Monte Carlo and prayer-bead Monte Carlo. 
We obtain extremely similar results concerning the tran- 
sit depth for the different wavelengths with the three tech- 
niques. The final results are listed in table [4] As reported by 
Carter al. (2008), there exists a degeneracy between the fit- 
ted orbital inclination i and a/i?*; whereas the ratio of radii 
k and the transit duration, T, are far more robust quantities. 
As a result of this robustness, we are able to use the fitted 
transit duration values as a test of whether the lightcurves 
appear physical or not. 

From our six fitted lightcurves, the duration of 5.8/im, 
epoch 1, cannot be used because only half the transit is 
fitted and so the fitted duration is dependent on priors. The 
other five lightcurves produce durations consistent with an 
average duration of T = 9508 ± 17 seconds with a Xred ~ 
2.6 suggesting an outlier. Removing the 3.6/im, epoch 2, 
measurement to leave us just the 4 preferred observations 
we find T = 9538±14 seconds with Xred = I- 1 - Since 3.6/im, 
epoch 2, produces an outlier duration and is also known to 
exhibit by far the strongest pixel phase effect out of all of 
the observed lightcurves, we give it zero weighting in the 
later spectral analysis. 

We find that our average durations are consistent with 
the duration we find when refitting the Brown et al. (2001) 
HST lightcurve of T = 9525l 14 seconds. Consequently, we 
conclude that our results support a solution consistent with 
the Brown et al. (2001) observations and thus we may be 
confident that the systematic corrections have been success- 
ful. 

4.9 Comparison with HD 189733b data 

Beaulieu et al. (2008) gives an accurate description of the 
method adopted to analyse the two IRAC channels at 3.6 
and 5.8 /jm in the case of the hot- Jupiter HD 189733b. The 
software BLUE (Alard, 2010), was used to fit the PSF, as 
several stars were in the field and could be used as calibra- 
tors. One of the capabilities of BLUE is to provide optimised 
centroid estimates, and provide an accurate modelling of the 



Water in HD 209458b 's atmosphere from 3.6—8 \im IR AC 'photometric observations in primary transit 



U.6W 





Beoulieu et ol., 2008 J 

Desert et al., 2009 J 















5000 10000 15000 

TIME (sec) 

Figure 7. The lower pannel shows the reprocessed HD 189733b 
data at 5.8 /im overplotted with the logarithmic correction from 
Desert et al. (2009) and the linear correction (Beaulieu et al. 
2008). The vertical line indicates 2000 sec. We provided in the 
text evidences to reject the first 2000 sec. This figure shows how 
the logarithmic correction is overcorrccting in the transit. 



PSF. In the case of HD 209458b one star only was present, 
so we had to adopt a different strategy, using the SExtractor 
programme (Bertin and Arnouts, 1996). The extracted light 
curves at 3.6 and 5.8 fim were corrected in a similar manner 
to that detailed in Beaulieu et al. (2008). In particular the 
3.6 ^m observations for the two planets show moderate or 
strong pixel-phase effects, that can be corrected for. 

5.8 lira observations ostensibly represents the greatest 
challenge for correcting systematic errors as the behaviour 
is somewhat less understood than the 8.0 fira ramp and the 
pixel phase effects at 3.6 fj,m and 4.5 /jm. This situation 
is exacerbated by the observation of a discontinuity in the 
photometry in two separate transit observations. At 5.8 /im 
there are no significant pixel-phase effects, but a linear drift 
with time, after the first 2800 s in the case of HD 209458b 
and 1800 s in the case of HD 189733b. For both planets we 
disregard the first 2800/1800 seconds and then simply apply 
a linear correction to the data after this point. 

This is the main discrepancy between the Beaulieu et 
al. (2008) reduction and the one adopted by Desert et al. 
(2009). We do not discuss here previous results or methods 
adopted by the same team, incorporated in Ehrenreich et al. 
(2007), in part because we have already explained the rea- 
sons of their inadequacy in the Beaulieu et al. (2008), but 
most importantly as clearly abandoned by the authors them- 
selves in the new version of the analysis of the same data 
provided by Desert et al. (2009). Desert et al. (2009) applied 
a logarithmic time-correlated detrending to this channel of 
the same form of 8.0 fim observations, but of opposite sign 
i.e. an 'anti-ramp'. In contrast, we applied a truncation of 
the first 2800 seconds followed by a linear time-trend de- 
correlation. 

The question as to which method is the correct one is 
naturally a topic of debate within the community but we 
believe we have produced here an in-depth analysis of each 
type of correction. In §3.2 we treated the 5.8/im data with 
both methods and compared the resultant lightcurves. Sev- 
eral tests suggest that the truncation + linear detrending 
produces a more physical transit signal and an improved 
overall fit. In this work, we acknowledge that there currently 
exists no widely-accepted physical explanation for the chan- 
nel 3 systematic effects and thus the preference between the 



Table 6. Comparison of values of transit depth for HD 189733b 
at 3.6 and 5.8 fim by Beaulieu et al, (2008), Ehrenreich et al., 
(2007), and Desert et al, (2009). 



IRAC Beaulieu 2008 Ehrenreich 2007 D esert 2009 

3.6 /mi 2.383 ± 0.014% 2.434 ± 0.026% 2.387 ± 0.0093% 
5.8/im 2.457 ±0.017% 2.375 ±0.04% 2.393 ± 0.016% 



linear and logarithmic model must be made primarily on the 
basis of the lightcurve information. On this basis, we cannot 
justify employing the logarithmic model over the method 
adopted here, given the range of evidences compiled. 

In our case, we find that adopting the logarithmic fit to 
our HD 209458 5.8/im data would underestimate the transit 
depth by 0.035 %, generating a systematic error of ~ 2.3<r. 
We also estimate that the incorrect use of logarithmic cor- 
rection leads to an under estimate of the transit depth of HD 
189733 by 0.047% , generating a systematic error of ~ 2.9a". 
This accounts for the discrepancy between the studies of 
Desert et al. 2009 compared to Beaulieu et al. 2008. How- 
ever, for case of 3.6 /im, both teams agree upon the correc- 
tive procedure, and so we should expect very similar results. 
Indeed, for the HD 189733 3.6/im photometry, the values and 
error bars estimated by Beaulieu (2008) and Desert (2009) 
are in excellent agreement, as shown in Table [(J 



5 DATA INTERPRETATION 

To interpret the data, we used the radiative transfer mod- 
els described in Tinetti et al. (2007a, b) and consider haze 
opacity, including its treatment in Griffith, Yelle and Mar ley 
(1998). 

Our analysis includes the effects of water, methane, car- 
bon dioxide, carbon monoxide, pressure-induced absorption 
of H2 — H2- We do not consider the presence of particu- 
lates, because there is no indication of particles large enough 
(~ 3/im ) to affect the planet's middle-infrared spectrum. 
The effects of water absorption are quantified with the BT2 
water line list (Barber et al., 2006), which characterises wa- 
ter absorption at the range of temperatures probed in HD 
209458b. Methane was simulated by using a combination of 
HITRAN 2004 (Rothman et al, 2005) and PNNL data-lists. 
Carbon monoxide absorption coefficients were estimated 
with HITEMP (Rothman et al. 1995) whilst for carbon diox- 
ide we employed a combination of HITEMP and CDSD-1000 
(Carbon Dioxide Spectroscopic Databank version for high 
temperature applications; Tashkun and Perevalov, 2008). 
The continuum was computed using H2 — H2 absorption 
data (Borysow et al., 2001). In fig. [8] we show the contribu- 
tion of the different molecules combined to water. 

The 3.6 /xm (and to a lesser degree the one at 8 /im) 
IRAC channel measurement can be affected by the presence 
of methane. By contrast, CO2 and CO may contribute in 
the passband at 4.5 fim. 

We find absorption by water alone can explain the spec- 
tral characteristics of the photometric measurements, which 
probe pressure levels from 1 to 0.001 bars (fig. I T0|) . The 
determined water abundance depends on the assumed tem- 
perature profile and planetary radius. We find that the data 
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Figure 8. Simulated middle Infrared spectra of the transiting Hot 
Jupiter HD 209458b in the wavelength range 3-10 fim. Water 
absorption is responsible for the main pattern of the spectra. The 
additional presence of methane, CO and CO2 are simulated in 
the blue, violet and green spectra respectively. 



Spectroscopic data are needed to further investigate the 
composition of this planetary atmosphere. 



6 CONCLUSION 

We have presented here IRAC photometry data record- 
ing the primary transit of HD 209458b in four infrared 
bands. We find that the systematics are very similar to those 
present in the data set obtained for the planet HD 189733b 
(Beaulieu et al., 2008), and therefore we adopted similar 
recipes to correct for them. We have performed Markov 
Chain Monte Carlo and prayer-bead Monte Carlo fits to 
the data obtaining almost identical results. Our observa- 
tions indicate the presence of water vapour in the atmo- 
sphere of HD 209458b, confirming previous detections of 
this molecule with different techniques/instruments. Inter- 
estingly, the thermal profiles derived for the day-side are 
compatible with this set of data probing essentially the plan- 
etary terminator. It is possible that additional molecules, 
such as methane, CO and/or CO2 are also present, but the 
lack of spectral resolution of our data have prevented these 
from being detected. Additional data in transmission at dif- 
ferent wavelengths and/or higher resolution will be required 
to gain information about these other molecules. 
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Figure 9. Modeled spectral absorptions of H2O in the atmo- 
sphere of HD 209458b for 1500, 2000 and 2500 K. 



can be interpreted, with a thermochemical equilibrium water 
abundance of 4.5 x 10" 4 (Liang et al., 2003;2004), assum- 
ing temperature profiles from Swain et al., 2009b. However, 
~ 1% difference in the estimate of the planetary radius, 
is compatible with water abundances 10 times smaller or 
larger, or with an overall change in the atmospheric tem- 
perature of about ~ 500K (see fig. [9]). Additional primary 
transit data at different wavelengths are needed to improve 
the constraint. 

While the contribution of other constituents is not nec- 
essary to interpret the measurements, mixing ratios of 10 -7 , 
10 -6 and 10 -4 of CO2, CH4 and CO, respectively, are al- 
lowed in our nominal model. 
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